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Abstract This article presents a guided introduction to a general class of interact- 
ing particle methods and explains throughout how such methods may be adapted 
to solve general classes of inference problems encountered in actuarial science and 
risk management. Along the way, the resulting specialized Monte Carlo solutions 
are discussed in the context of how they compliment alternative approaches adopted 
in risk management, including closed from bounds and asymptotic results for func- 
tionals of tails of risk processes. 

The development of the article starts from the premise that whilst interacting particle 
methods are increasingly used to sample from complex and high-dimensional dis- 
tributions, they have yet to be generally adopted in inferential problems in risk and 
insurance. Therefore, we introduce in a principled fashion the general framework of 
interacting particle methods, which goes well beyond the standard particle filtering 
framework and Sequential Monte Carlo frameworks to instead focus on particu- 
lar classes of interacting particle genetic type algorithms. These stochastic particle 
integration techniques can be interpreted as a universal acceptance-rejection sequen- 
tial particle sampler equipped with adaptive and interacting recycling mechanisms 
which we reinterpret under a Feynman-Kac particle integration framework. These 
functional models are natural mathematical extensions of the traditional change of 
probability measures, common in designing importance samplers. 
Practically, the particles evolve randomly around the space independently and to 
each particle is associated a positive potential function. Periodically, particles with 
high potentials duplicate at the expense of low potential particle which die. This 
natural genetic type selection scheme appears in numerous applications in applied 
probability, physics, Bayesian statistics, signal processing, biology, and information 
engineering. It is the intention of this paper to introduce them to risk modeling. 
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1 Introduction to Stochastic Particle Integration 

The intention of this paper is to introduce a class of stochastic particle based in- 
tegration techniques to a broad community, with a focus on risk and insurance 
practitioners, who we believe will benefit from the development of such inferen- 
tial methods in several problems encountered in these domains. A key motivation 
for this endeavor is due to the fact that over the last two decades, stochastic particle 
integration models have been extensively and routinely used in engineering, statis- 
tics and physics under sometimes different names, such as: particle filters, bootstrap 
or genetic filters, population Monte Carlo methods, sequential Monte Carlo mod- 
els, genetic search models, branching and multi-level splitting particle rare event 
simulations, condensation models, go-with-the winner, spawning models, walkers 
population reconfigurations, pruning-enrichment strategies, quantum and diffusion 
Monte Carlo, rejuvenation models, and many others. They have however not yet 
been routinely applied to develop solutions in important financial domains such as 
those we discuss in this tutorial type review, in which we discuss and illuminate 
areas they will benefit the fields of risk and actuarial science. To contain the scope 
of such an endeavor we focus the application domain discussions on the widely 
used class of problems based around the consideration of single risk processes un- 
der the Loss Distributional Approach, see detailed discussions in HI, |l2l,|[3l and in 
the dynamic setting ID and ||5l . 

We begin with the introduction of the fundamental background for interacting 
particle systems highlighting key papers in their developments through a range of 
different science disciplines, before introducing aspects of these stochastic methods 
to risk and insurance. It is important that practitioners when first encountering such 
Feynman-Kac interacting particle methods are aware that they encompass a far more 
general gamut of stochastic integration and optimization methodologies then the 
most well known sub-class of such methods knows as particle filters, which are 
typically utilized in inference under latent process state space models in engineering 
and statistics. It is the intention of this article to explain the key papers and ideas in 
a far more general framework which is much more encompassing than the special 
subset of particle filter based algorithms. It is for this reason that me mention that 
the stochastic methods discussed in this paper are applicable to a significantly larger 
subset of problems than the standard particle filter approach that may be thought of 
when first encountering Sequential Monte Carlo. 

Then we proceed through a selection of key features of their development, focus- 
ing on a sub-class of such methods of relevance to the application domain explored 
in this manuscript, risk and insurance. The origins of stochastic particle simulation 
certainly starts with the seminal paper of N. Metropolis and S. Ulam 16]. As ex- 
plained by these two physicists in the introduction of their pioneering article, the 
Monte Carlo method is, "essentially, a statistical approach to the study of differ- 
ential equations, or more generally, of integro-differential equations that occur in 
various branches of the natural sciences". The links between genetic type particle 
Monte Carlo models and quadratic type parabolic integro-differential equations has 
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been developed in the beginning of 2000' in the series of articles on continuous time 
models Qll. 

The earlier works on heuristic type genetic particle schemes seem to have started 
in Los Alamos National Labs with works of M.N. Rosenbluth and A.W. Rosen- 
bluth 191, and T.E. Harris and H. Kahn ifTOl . We also quote the work on artificial 
life of Nils Aall Barricelli lfmfT2l . In all of these works, the genetic Monte Carlo 
scheme is always presented as a natural heuristic resampling type algorithm to gen- 
erate random population models, to sample molecular conformations, or to estimate 
high energy particle distributions, without a single convergence estimate to ensure 
the performance, nor the robustness of the Monte Carlo sampler. 

The mathematical foundations, and the performance analysis of all of these dis- 
crete generation particle models are rather recent. The first rigorous study in this 
field seems to be the article jlTl published in 1996 on the applications of particle 
methods to nonlinear estimation problems. This article provides the first proof of 
the unbiased property of particle likelihood approximation models (lemma 3 page 
12); and adaptive resampling criteria w.r.t. the weight dispersions (see remark 1 on 
page p.4). We also quote the first articles presenting heuristic type particle filters 
lfT4l[T5]| . and a series of earlier research reports lfT6l[T7l[T8l[T9l . 

For an in depth description of the origins of particle methods, and their appli- 
cations we refer to the following studies ll20l |2T1 . These articles also contain new 
stochastic models and methods including look-ahead type strategies (section 4.2.2), 
reducing the variance using conditional explorations w.r.t. the observation sequences 
(example 3 p. 40), local errors transport models (see the proof of theorem 1 on page 
1 1), mean field models w.r.t. the occupation measures of random trees (section 3.2). 

A more detailed review of particle models in discrete and continuous time can 
be found in 12211231 . In the research monograph the reader will find a detailed dis- 
cussion on particle models and methods including acceptance-rejection with recy- 
cling particle strategies, interacting Kalman filters a.k.a. Rao-Blackwellized par- 
ticle filters (section 2.6, and section 12.6.7), look-ahead type strategies (section 
12.6.6), genealogical tree models and branching strategies (section 11), and inter- 
acting Metropolis-Hasting models (chapter 5). 

The practitioner will find in the research books 1231 1241 a source of useful con- 
vergence estimates as well as a detailed list of concrete examples of particle ap- 
proximations for real models, including restricted Markov chain simulations, ran- 
dom motions in absorbing media, spectral analysis of Schrodinger operators and 
Feynman-Kac semigroups, rare event analysis, sensitivity measure approximations, 
financial pricing numerical methods, parameter estimation in HMM models, island 
particle models, interacting MCMC models, statistical machine learning, Bayesian 
inference, Dirichlet boundary problems, nonlinear filtering problems, interacting 
Kalman-Bucy filters, directed polymer simulations, stochastic optimization, and in- 
teracting Metropolis type algorithms. 

There is an extensive number of texts on particle simulation and sequential 
Monte Carlo samplers, many of them contain much practically oriented discus- 
sions including Bayesian inference, nonlinear filtering and optimization, as well 
as optimal control problems. For a further discussion on the origins and the ap- 
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plications of these stochastic models, we refer the reader to the following texts 
|26l 1221 |28l 123 [3i ini [331 [31 13^ references therein. 

Particle methods are yet to be routinely or widely introduced to areas of risk and 
insurance modeling. The initial examples that have been developed are detailed in 
ll36l . where a special sub-class of such methods was developed for an important set 
of risk management problems. It is the intention of this paper to highlight aspects of 
this class of problems and the stochastic particle solutions. 

Before, proceeding with the introduction of stochastic particle methods, we first 
provide an overview of the domain of application to be considered in risk and insur- 
ance and a mathematical discussion that will motivate the relevance of such stochas- 
tic particle integration methods. 



2 Motivation for Stochastic Particle Solutions: 

How Such Methods May Complement Risk Process Asymptotics 

Here we provide motivation to explain how and why risk management and actuarial 
sciences can benefit from the development of interacting particle system inferential 
solutions to an important subset of generic problems faced by practitioners in these 
domains. In particular we focus on a an aspect of single risk loss processes described 
under a Loss Distributional Approach (hereafter LDA) framework, see discussion 
in II37I . 0, Q and the books 1381 and [?] for the background on such modeling 
approaches in risk. For basic discussions on how such problems relate to a large 
class of non-life insurance problems see examples in 121 . 



2.1 The Loss Distributional Approach and Risk Management: 
a tale of light to heavy tails 

In this section we fist motivate and introduce the context of LDA modeling in risk 
and insurance. Then we present three key challenges associated with working with 
such models faced by risk and insurance practitioners, thereby effectively detailing 
important inference challenges faced by such practitioners. Next, we provide a brief 
specifically selected survey of closed form analytic results known in the actuarial 
and risk literature for sub-classes of such LDA models as the Single Loss Approx- 
imations (hereafter SLA). We detail the closed form solution for the light-tailed 
severity distribution case and then explain how such approaches can not be obtained 
in such a form in the heavy-tailed sub-exponential risk process settings, often of in- 
terest in the domain of risk and insurance. As a result, we briefly present the results 
recently developed in actuarial literature for the heavy tailed case corresponding to 
the first order and second order asymptotic approximations, see comprehensive dis- 
cussions in a general context in l39l , l40l and the books, 11411 and the forthcoming 

ia. 
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We conclude this section by observing that according to regulatory standards 
and indeed good risk management practice such approximations are often required 
to be accompanied with numerical and statistical solutions which can more readily 
take into account model uncertainty, parameter uncertainty and the fact that such 
approximations are inherently asymptotic in nature, and may be inaccurate outside 
of the neighborhood of infinity. In this regard we summarize a class of interacting 
particle solutions recently developed to address such estimations which complement 
such closed form asymptotic results. 

Consider the widely utilized insurance model known as a single risk Loss Distri- 
butional Approach model. This represents the standard under the Basel Willi capital 
accords 143)1 and involves an annual loss in a risk cell (business line/event type) 
modeled as a compound distributed random variable, 

zl^^^l^x^Ht), (1) 

.s = l 

for f = 1 , 2, . . . , r discrete time (in annual units) and index j identifies the risk cell. 
Furthermore, the annual number of losses is denoted by A^/^' which is a random 
variable distributed according to a frequency counting distribution P'^'(-), typically 
Poisson, Binomial or Negative Binomial and the severities (losses) in year t are 
represented by random variables ^'^'(f), i > 1, distributed according to a severity 
distribution F^j^-). 

In constructing this model we assume that all losses are i.i.d. with Xp^ (f ) - 
Fx{x) and that the severity distribution is continuous with no atoms in the support 
[0,oo). As a consequence, linear combinations (aggregation) of losses in a given 
year, denoted by 

Sit,n)^txPit)^Fs{x) 
have the following analytic representation: 

Fsix)^iF*F*---F){x)^ [ F^"-^'>*{x-y)dF{x). 

i[0,c») 

We also observe that due to a result in ||44| if F{x) has no atoms in [0,0°) then the 
n-fold convolution of such severity distributions will also admit no atoms in the 
support [0,oo). The implications of this for such Interacting Particle based numeri- 
cal procedures (IS, SMC, MCMC) is that it ensures numerical techniques are well 
defined for such models when considering ratios of densities on the support [0,oo). 
In addition we note that continuity and boundedness of severity distribution Fx{x) 
is preserved under «-fold convolution if Fx{x) admits a density j^Fx{x) then so 
does Fs{x). For most models such analytic representations of the combined loss dis- 
tribution are non closed form, with the exception of special sub-families of infinitely 
divisible severity distribution models, see B31 . 
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It is important in practice to consider carefully the single risk processes in which 
business managers believe will produce infrequent losses with very high conse- 
quence. Modeling such risk processes typically requires sub-exponential severity 
distributions. If one considers losses , . . . ,X„,. . . as independent positive random 
variables with distribution F (x) = P (Xj. < x) , e { 1 , 2 ,...,«,...} . Then the class 
of sub-exponential distributions (F{x) G J^) satisfy the limits 

limi^=. (2) 

if and only if 

\-F^-Hx) 

hm — — Y = 2. (3) 
The sub-exponential family of distributions ^ defines a class of heavy-tailed sever- 



ity models in which B6I demonstrated the necessary and sufficient condition for 
membership being, for a severity distribution F g ^ if and only if the tail distribu- 
tion F(x) = 1 — Fix) satisfies 

rF(x-y) , , 
lim / -^^^F{j)dy=\. 

Alternatively one may characterize the family of distributions F E -J^ hy those that 
satisfy asymptotically the tail ratio 

lim^^|-^ = l,Vye[0,oo). (4) 

Severity models F € ,^ are of interest for severity distributions in high consequence 
loss modeling since they include models with infinite mean loss and infinite vari- 
ance. In addition, the class ^ includes all severity models in which the tail distri- 
bution under the log transformed r.v., F (log(x)), is a slowly varying function of x at 
infinity. 

To further understand LDA modeling with sub-exponential severity models we 
recall the notion of asymptotic equivalence in which a probability distribution func- 
tion F{x) is asymptotically equivalent to another probability distribution function 
G{x), denoted by F(x) ~ G(x) as a: — >■ 0° if it holds that, Ve > 0, 3xq such that Vjc > xq 
the following is true 

G{x) 

Furthermore, we say that a probability distribution function is max-sum-equivalent, 
denoted hy F G, when the convolution of the tail distribution of two random 
variables is distributed according to the sum of the two tail distributions asymptoti- 
cally, 

1 - iF*G){x) = (F*G)(x) - F{x) +G{x), x^oo. 



-1 



< e as X ^ 00 (5) 
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Then for the class of heavy tailed sub-exponential LDA models we have that a prob- 
abihty distribution function F will belong to the sub-exponential class ^ifF^Mp, 
i.e. it is max-sum-equivalent with itself and that the class ^ is closed under convo- 
lutions. The implications of this for the LDA models are clear when one observes 
that sub-exponential LDA models are compound process random sums comprised 
of an infinite mixture of convolved distributions, 

G(x) = £A„r'*W, (6) 

n=0 

for a suitable series {A„}, (e.g. convergent sequence satisfying Kolmogorov three 
series theorem). 

Proceeding in this section we will consider the compound Poisson distribution 
case in which the sequence A„ = e^^^ and one can show the practically relevant 
asymptotic equivalence between the severity distribution F and the annual loss dis- 
tribution G such that if F e then Ge ^ and 

lim£|4=A. 

These properties of sub-exponential LDA models regarding asymptotic equiva- 
lences make stochastic quantile and tail expectation approximations tractable for 
estimation. Only is special families of infinitely divisible severity models can closed 
form annual loss distributions be obtained, see discussions in BtI and 1451 . 

In general based on these properties we can obtain asymptotic approximations 
to the annual loss distribution tails which typically fall under one of the following 
classifications: 

• "First-Order" and "Second-Order" Single Loss Approximations: recently dis- 
cussed in BSl . BOl . P9l and references therein. 

• "Higher-Order" Single Loss Approximations: see discussions in ll50l and recent 
summaries in 1391 and references therein. 

• Extreme Value Theory (EVT) Single Loss Approximations (Penultimate Ap- 
proximations): the EVT based asymptotic estimators for linear normalized and 
power normalized extreme value domains of attraction were recently discussed 
in 1491. 

• Doubly Infinitely Divisible Tail Asymptotics given a-stable severity models dis- 
cussed in 145] and g?! 

We now briefly detail the first and second order asymptotics known for light 
and heavy tailed severity distributions in LDA models, before then explaining how 
stochastic particle methods can be utilized to complement such closed form expres- 
sions and the role we believe will be played in the future by such stochastic algo- 
rithms in complementing these results. 
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2.1.1 The light tale of Ught tails 

At this stage it is very informative to understand the asymptotic results known for 
light tailed models as this will inform the results obtained in the heavy tailed ex- 
pansions. An elegant summary of such a result was provided by ISTl where they 
consider frequency distributions p„ = Pr(A^ = «) satisfying 

Pn ^ w"n^C(n), n — > oo. 

for some we (0, 1), y £ E and a function C(«) slowly varying at oo. If fc > exists 
such that the Laplace transform of the severity 

Lx{s) ^^[F{x)] = / exp(-M)t/F(x), G M, 





matches the radius of convergence of the generating function of the frequency dis- 
tribution, 

w"' =Lx{-k) 

with —L'^{—k) < oo, then the following asymptotic equivalence for the compound 
process tail distribution is satisfied. 

This light tailed asymptotic result demonstrates that the behavior of the com- 
pound loss distribution tail is determined by either the frequency or the severity 
depending on which has the heavier tail. In addition it is clear that the Poisson dis- 
tribution tail is too light for this result to be valid since the radius of convergence of 
generating function is infinite. There are therefore alternative expansions developed 
for compound Poisson risk processes such as the Saddle point approximation, see 
[]. If the severity distribution is bounded, then as x °o 

W - , ,..T!!!""^'^^l/2 [exp(-A(l -Lx(Jc)) -exp(-A)] , 

where K is the solution of —XL'^{k) ~ x. 

So how do these results relate and motivate the context we are considering in 
sub-exponential LDA models? 

Quite simply, in the sub-exponential heavy tailed setting the Laplace transform 
does not exist and hence these results do not apply. Examples of such models for 
which this is true include severity distributions with power law tail decay {Pareto, 
Burr, log gamma, Cauchy, a-Stable, tempered stable and t-distribution). 
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2.1.2 The heavier tale of heavy tails 

In this subsection we detail briefly the asymptotic fist and second order tail results 
for the LDA models when sub-exponential severity distributions are considered. The 
sub-exponential LDA first order tail asymptotics involve obtaining the closed form 
approximate expression for Fzf^{x), see details in ||52l . ||49l . To proceed, consider 
the annual loss distribution G{z) = Fz{z) under LDA formulation, given by, 

G{z) = Fz{z) = ^ [Z < z\N = n]Pr[N = n] = PnF^"^*{z), 

with the severity distribution satisfying Fx{z) G Furthermore, assuming that for 
some e > 0, 

£(1 

Then the right tail of the annual loss distribution Fz{z) for the annual loss random 
variable Z, is approximated according to a SLA given by, 

F^{x)^E[N]F^{x){l+o{l)) asjc^oo, 

or equivalently the tail of the annual loss distribution is sub-exponential in behavior 
with asymptotic equivalence, 

Fk{x)^E[N]Fx{x), 



To understand the basic result of the first order tail asymptotic Fz^ (x) consider 
two steps: 

1. Obtain an upper bound on the asymptotic ratio of Fz„{^) ™d severity F{x) for 
all n £ J. Typically one can apply Kesten's Bound which states that for subexpo- 
nential severity distributions F there exists a constant K = ^(e) < °° for e > 
s.t. V« > 2 the following bound holds dlSSl ) 

-^<K{l+8y\x>Q. 
F{x) 

2. Then simply utilize the Kesten bound to motivate the application of dominated 
convergence theorem to interchange the order of summation and limit and recall 
characterization of heavy-tailed sub-exponential severity models. 



,. F*^{x) ^ . ,. ,. F*"(x) 
lim ^= = 2, implies lim ^= 

x^oo f(^x) ^■^"o F{x) 



This process gives Fzf^{x) ^ E[N]F{x) since: 
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lim %^ = lim £ p,,^ = £ np„ = ]E[A^], 



As discussed in |40j, and the papers therein, the second order asymptotic results 
can be developed in a wide class of risk models by considering the following further 
assumptions. 

Assumption 1: F is zero at the origin (x = 0) and satisfies that both the tail 
distribution F and density f are subexponential. 

Assumption 2: The frequency distribution N ^ FN{n) is such that its probability 
generating function given by 



p^(v)=E[v^] = £Pr(A^ = n)v", 

n=0 



is analytic at v = 1 



Examples of severity models that satisfy such assumptions include: Log-Normal, 
Weibull (heavy tailed), Benktander Type I and Type II, Inverse Gaussian, a-Stable, 
Halphen Family, Normal-Inverse-Gaussian and other members of Generalized Hy- 
pergeometric family. 

When distributions satisfy Assumption 1 and Assumption 2 then two situations 
may arise, those in which the loss random variable distribution is finite mean and 
the alternative case in which an infinite mean loss model is considered. If the loss 
r.v. has finite mean {E[X] < °o) then the following result can be derived, see 1541 
and ll55ll for details. 

^^FzW EWFW^^^ (7) 
f[x) 

Alternatively, if the loss r.v. is infinite but the severity density satisfies / e 
for I < P <°° then: 

,. Fz(x)-E\N]F(x) , , 

f{x)jQF{s)ds 

-T^f 1 — 1/6) 

with ci = 1 and = (1 ^ 1^) 2r\i-2/p) ^'^^ ^ ^ (1'°°)- We recall that a measurable 
function / : (0, °o) — > (0, 0°) is called regularly varying (at infinity) with index /3 and 
denoted by / e if Vm > 0, 3)3 e K where 

fiux) r. 
hm = m'^ . 

x^»o f(x) 

In the following subsection we clearly detail how and why such asymptotic re- 
sults are utilized for inference in risk and insurance, before highlighting the im- 
portant potential role stochastic particle methods will play in complementing these 
results. 
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2.2 Inferential Challenges for Risk and Insurance: 
Asymptotics and the Role for Stochastic Particle Integration 

The asymptotic approximation methods just surveyed were developed in the actu- 
arial literature to tackle the serious statistical and computational challenges posed 
by estimation of tail quantiles and expectations for heavy-tailed LDA models. The 
continued interest in such asymptotic results primarily stems from the fact that such 
closed form expressions bypass the serious computational challenges for estimation 
of risk measures for such heavy-tailed annual loss distributions under traditional 
integration methods, Fourier methods, recursions (Panjer) or basic Monte Carlo ap- 
proaches. However, they do have associated issues, see discussions in ||56l . 

The properties of such asymptotic single loss approximation estimates are still an 
active subject of study with regard to aspects such as explicit approximation errors, 
unbiased quantile function estimation, asymptotic rates of convergence, sensitivity 
to parameter estimation and model misspecification. It is in the understanding of 
these increasingly important practical features that we believe stochastic particle 
integration methods will complement the asymptotic results, which are generally 
un-attainable under such additional considerations. 

Before introducing in depth some key results in stochastic particle methods, we 
will first complete the risk and insurance motivation by explaining the key result 
one obtains from a risk management perspective as a consequence of these first and 
second order asymtptotics. We will also demonstrate exactly how it is often utilized 
to make risk management based decisions by tying it back to the calculation of risk 
measures in LDA models. 

Based on the results obtained for the second order asymptotic in the heavy tailed 
LDA models, one can show that if the severity distribution F satisfies Assumption 1 
and Assumption 2 with a finite mean, and the hazard rate h{x) ^ ^ T is of regular 
variation h € RV_^ for j3 > 0, then as a 1 one has for the inverse of the annual 
loss distribution the result 



From this result it is then possible to consider asymptotic approximations of key risk 
management quantities known as risk measures which are used in the allocation of 
capital and reserving in all financial institutions and stipulated as standards under 




(8) 



where a = 1 - (1 - a)/¥.[N] and 
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regulatory accords in both Basel II/III and Solvency II. Examples of such tail func- 
tionals include the calculation of Value-at-Risk (VaR), Expected Shortfall (ES) and 
Spectral Risk Measures as detailed below in both their definitions and the resulting 
simple asymptotic approximations one may consider. 

These asymptotic expansions allow one to obtain estimates of common risk mea- 
sures, see 1571 and JH], such as Value-at-Risk (VaR) for a level a e (0, 1 ), given by 
the quantile of the annual loss distribution, 

VaRz (a) = F^{a) = inf {z e M : Fziz) > a} 



where F^{-) is the generalized inverse, see ||59l . The Expected Shortfall (ES), see 
ll60ll . for a level a £ (0, 1) is given by the tail expectation of the annual loss distri- 
bution according to 

ESz(a) =E[Z|Z> VaRz(a)] = / YaRzis)ds 

, , . (10) 
a , f 1 - a \ a 

~ F^ 1 r-r VaRz (a) , 

a-l \ E[N] J a- 1 ^ 

and the Spectral Risk Measure ( SRM) for a weight function : [0,1] i— > R given by 
SRMz(0)= ^{s)WaRz{s)ds 



^J(f{a,^i)F^ ^l--^j^.jf(a,(^i)VaRzia), 

with Vf e (l,oo) a function 0i(l - 1/f) < Kt'^^l^+^-'^ for some AT > and e > 
where ^ 

je'iaji)^ s^/P--(l)iil-l/s)ds. 



2,2.1 The Role for Stochastic Particle Methods 



Though the asymptotic results presented are elegant and efficient to evaluate, they do 
warrant careful consideration in their application. In this section we explain what we 
mean by this statement and then utilize this to motivate the use of stochastic particle 
methods. It is important for practitioners to understand that the properties of such 
SLA estimates is still a subject of study, this includes for example an understanding 
of (approximation error, unbiased quantile function estimation, asymptotic rates of 
convergence, sensitivity to parameter estimation, model misspecification etc.) all of 
which can have a non-trivial influence on the resulting asymptotics and therefore 
risk measure estimates, as discussed recently in f6l\ . 
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In practice it may often be the case that one requires calculation of VaR, ES and 
Spectral Risk Measures at levels which do not satisfy such asymptotic properties, 
rendering such approximations inaccurate. In addition, though not yet a regulatory 
requirement, it is always good practice to consider the uncertainty associated with 
the estimation of the tail functionals and quantiles, through say confidence intervals 
and this is non-trivial to obtain under such asymptotic expansion results. Thirdly, as 
discussed in 1391 and BTI the asymptotic rates of convergence of such approxima- 
tions are still only known in a little-oh Landau sense and therefore do not inform or 
guide the applicability of such results. There is also a significant complication with 
such asymptotic results that arises in the application of such results for models that 
in practice do not admit closed form representations of the quantile function of the 
severity distribution in the LDA model. Finally, there is also a significant interest in 
diversification benefits that may be gained through the modeling of tail dependence 
features in the multi-variate risk process setting. Whilst these asymptotic results are 
extendable to the multi-variate case of multiple risk processes, the addition of even 
parametric tail dependence through a copula renders the derivation of such results 
highly challenging and an active research area at present, see for example recent 
results in 1621 . 

It is in these four key elements that we argue stochastic particle based numerical 
solutions to such inference on risk measures and tail functionals can be of direct 
utility to complement such asymptotic results. However, as all practitioners will 
know, the naive implementation of standard Monte Carlo and stochastic integration 
approaches to such problems will produce often poor results even for a considerable 
computational budget, see discussions in f63l/ . There is therefore a computational 
challenge for estimation of risk measures for such heavy-tailed annual loss distribu- 
tions that we argue can be addressed by stochastic particle methods. This concludes 
the discussion on motivations for how such methods can play an increasingly more 
important role in risk and insurance and we now present a detailed exposition of 
a few important stochastic particle methods that will be a solid starting point for 
practitioners. 



3 Selected Topics in Stochastic Integration Methods 

In this section we will introduce a variety of stochastic integration methods, present- 
ing them formally from a mathematical perspective and making clear the properties 
of such methods. This will provide practitioners with an understanding of the key 
properties of these methods and the relevant references to consider in applying such 
approaches to tackling risk and insurance problems. Note, in this section the notation 
adopted is utilized to reflect that which is considered in the statistics and probability 
literature where much of the formal study of these methods has taken place. We first 
introduce examples of problem domains where each approach has been considered 
in the risk and insurance literature, to tackle particular inference problems for spe- 
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cific models, before formally detailing the general mathematical understanding of 
these methods. 



3.1 Standard Monte Carlo Techniques for Risk and Insurance 

Let Z be a li-dimensional random variable and A C M'' some measurable subset. 
Suppose we want to compute the quantity F{X E A) := Fx{A). We further assume 
that it is straightforward to generate a sequence {X')i<j<N of independent copies of 
the random variable X. In this situation, the traditional Monte Carlo approximation 
of the distribution Fx is given by the empirical measures 

^x = ^ L Sx.^Nt^Vx 

" 1<(<A' 

More precisely, the convergence can be understood as the weak convergence of 
empirical measures, in the sense that the following convergence holds 

P^(/):=//(x)IPf(^x) = l £ fir)^^^^Fxif)= f fix)Fx{dx)=E{f{X)) 

•' \<i<N •> 

almost surely, for any bounded measurable function / on . Using indicator func- 
tions of cells in M'', the shape of the measure Fx can be obtained by plotting the 
histograms of the samples X' in every dimensions. By the strong law of large num- 
bers, the above convergence is also met for integrable functions w.rt. the measure 

For indicator functions / = 1a, sometimes we make a slight abuse of notation and 
we set P^(A) and Px(A) instead of P^^(Ia) and Px(1a)- From the above discussion, 
we already have that 

^(^):=77 E ^ivt~Px(A)=E(U(X)). 

" l<i<A' 

The following properties are readily checked 

E(P^^(A)) = Px(A) and Var (Pf (A)) = i Px(A) (1 - Pz(A)) . 

In addition, an A^- approximation of the conditional distribution of X w.rt. the event 
{X E A} is given by 

^ lAix)F^idx)~>N^„—^lAix)Fxidx)=FiXedx\XeA). (12) 



P^^(A) ' ' ' Px(A) 



The l.h.s. terms in the above display is well defined as soon as F^{A) > 0. For 
rare event probabilities Px(A), say of order 10^^, the practical implementation of 
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this Monte Carlo algorithm meets the difficulty that we need too many samples to 
estimate Px{A) using the proportion of success of such an event occurring only once 
per millions of attempts. 

We illustrate this on a standard model in risk and insurance based on the Poisson- 
Log Normal LDA model of a single risk process. This example though simple is 
both widely utilized in practice and also illustrative of the complementary role of 
the asymptotic approximations and the role Monte Carlo plays, since this specific 
model admits a closed form expression for the survival quantile of the annual loss 
under the first order asymptotic. 

Example 1 (Single Risk LDA Poisson-Log-Normal Family). Consider the heavy- 
tailed severity model, selected to model the sequence of i.i.d. losses in each year r, 
denoted {^i(0}(=i Ar,' ™d chosen to be a Log-Normal distribution X, ~ LN{^,a) 
where the two parameters in this model correspond to parameterizing the shape 
of the distribution for the severity a and the log-scale of the distribution /i. The 
survival and quantile functions of the severity are given by 

fx{x;ii,a) = — e 2ci2 , X > 0; jLi e M (7 > 
xwlita^ 

F{x\jx,a) = 1 - F{x) = ^^== exp [^-^^ Oog(«) - m')) du 
2 • 2"' ,^>0;jueIR(T>0 




e(p)=exp(jU + (J4' '(p)), 0<p< 1. 

"ore the closed form 5 
in this case under a fi 
according too Equation [T3] 



Therefore the closed form SLA for the VaR risk measure at level a would be pre- 

■^N - 



sented in this case under a first order approximation for the annual loss Z = Y^=\^i 



VaRa [Z] = exp 



1-a 



(13) 



To compare this first order asymptotic result to the crude Monte Carlo approach 
(for which one can generate uncertainty measures such as confidence intervals in 
the point estimator) it is first required to detail how to simulate such an annual loss 
process. In this simple example the simulation of a loss process from a Log-Normal 
severity distribution can be achieved via a transformation of a Gaussian random 
variate, which itself is generated typically via a transformation of two uniform ran- 
dom variates in a Box-Muller approach, see details in I64II . The basic Monte Carlo 
simulation of the annual loss process for T-years, from a Poisson-Log-Normal LDA 
model can be achieved via a transformation of standard random variates as follows: 

Algorithm 1: Poisson-Log Normal LDA model via Standard Monte Carlo 

1. Generate vector of realized annual loss counts Ni-j = {Ni ,N2, ■ ■ ■ ,Nt} by draw- 
ing from a Poisson distribution with rate A, Nt ^ Po{X). This is undertaken 
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for each random variate realization A', = n, for each year of simulation t G 
{1,2, . . . , r} via one of several possible algorithms, such as: 

a. Set L = exp(— A ) , ^ = and p = I- 

b. While p > L: set i = ,v+ 1 then generate U ^ U[0, 1] and set p = pU. 

c. Set n, = s 

2. For each t, generate «r i-i-d reahzations of the loss X,(f ) ~ LA^(/i, a) via transfor- 
mation according to: 

a. Generate two independent uniform random numbers Ui, U2 uniformly dis- 
tributed on (0, 1] 

b. Apply the Box-Muller transformation 

Y = IX + a (^^/T^\MU2) cos {InUi)^ {11, a^) . (14) 

c. Apply the transformation from Normal to Log-Normal 

X = exp(y)-LA^(jU,(7). (15) 

3. Set the realized annual loss in each year t to Z(f ) = L/^i ^KO- 

To complete this example, we illustrate the basic Monte Carlo solution for the VaR 
for a range of quantile levels of the annual loss distribution along with the measures 
confidence intervals in the point estimators, compared to the first order asymptotic 
result. The quantiles a G {0.70,0.75,0.80,0.85,0.9,0.95,0.99,0.995,0.9995} are 
considered where the 99.5% and 99.95% quantile levels do in fact correspond to 
regulatory standards of reporting in Basel II/III. 

This example illustrates clearly the motivation for consideration of more devel- 
oped methods, since one can see that even in this relatively simple example, de- 
pending on the values of the parameters in the LDA risk model, the asymptotic VaR 
approximation may or may not be accurate at quantile levels of interest to risk man- 
agement. Clearly in the top sub-plot the accuracy of a capital estimate formed from 
the SLA will be poor where as the accuracy from the second subplot with the differ- 
ent risk profile parameters is fine for quantiles beyond the 80-th percentile. However, 
in general since the rate of convergence is still an active topic of research for such 
approximations, the only way to ensure accuracy of such methods for a given set of 
estimated / specified parameters in practice is to complement these approximations 
with a numerical solution or comparison. In general such an approach will require 
more sophisticated numerical procedures for more challenging risk process settings. 
In this example we utilized 1,000,000 annual years and the Mote Carlo accuracy was 
sufficient, note as X decreases, as will be the case for high consequence loss events 
that one considers modeling with sub-exponential models, this number of samples 
will need to significantly increase. 
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— VaR Monte Carlo 

VaR Monte Carlo 95% Upper CI 




Quanllle Level for Value-at-Risk - Capital Level 




Quanlile Level for Value-al-Rlsk - Capital Level 



Fig. 1 Annual Loss VaR Capital Estimate versus quantile level for Poisson-Log Normal LDA Risk 
Process. Top Plot: Severity distribution /i=l,CT = 0.5,A = 3. Bottom Plot: Severity distribution 
jU = 1,CT = 5, A = 3. 

3.2 Importance Sampling Techniques for Risk and Insurance 

One could argue that the second most widely utilized class of stochastic integration 
methods utilized in risk and insurance settings would have to be the Importance 
Sampling family, see for example reviews with control variates in market risk set- 
tings in 1651 . In the insurance setting one can see for example 1361 . 

To understand these class of methods as another way to design a feasible algo- 
rithm is to sample using another random variable for which the occurrence prob- 
ability of the desired event P(y e A) := P}'(A) is closer to 1. This well known 
importance sampling strategy often gives efficient results for judicious choices of 
twisted measures Py. Nevertheless, in some practical situations, it is impossible to 
find a judicious Py that achieves a given efficiency. Furthermore, this importance 
sampling technique is intrusive, in the sense that it requires to change the reference 
statistical or physical model into a twisted sampling rule. 

To be more precise, sampling independent copies (F')i<,<Af with the same 
dominating probability measure Py ^ Px, the traditional Monte Carlo approxima- 
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tion is now given by 




= Px(A) 



The following properties are readily checked 
and 

It is easily check that 

Md^) - ^ F.(^x) ^ Var (f^ (l, ^)) ^ 

In other words, the optimal twisted measure Py is given by the unknown conditional 
distribution of X w.rt. the event {X G A}. In practice, we try to find a judicious 
choice of twisted measure that are easy to sample, with a probability mass distribu- 
tion that resembles as much as possible to the desired conditional distribution. 

Another traditional idea is to use the occupation measure of a judiciously cho- 
sen Markov Chain Monte Carlo (abbreviate MCMC) sampler with prescribed target 
measure 

7]{dx):=P{X edx\X eA). 

Of course, the first candidate is to take a sequence of independent copies of 
random variables with common distribution 77 . Several exact sampUng techniques 
can be used, including the inversion of the repartition function, change of variables 
principles, the coupling from the past, and acceptance-rejection techniques. For in- 
stance, the Monte Carlo approximation presented in (fTSl i is clearly based on this uni- 
versal and traditional acceptance-rejection sampling technique. A random sample Xt 
with distribution Px is accepted whenever it enters in the desired subset A. In this in- 
terpretation, we need to sample independent copies of X to obtain := x P^(A) 
independent samples with common law r\. However, for probabilities Px(A) of or- 
der 10^^, this method requires millions of samples. 



3.3 Markov chain Monte Carlo for Risk and Insurance 

After the class of standard Monte Carlo and the Importance SampUng techniques 
one would naturally consider the next most widely developed class of methods as 
the Markov Chain Monte Carlo methods. These have found use in insurance appli- 
cations in non-life reserving models for example in Chain Ladder models I166I . Il67ll . 
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||68ll and Paid Incurred Claims models ||69l and llTOl . in Operational Risk models in 
El, HI and in credit risk modeling for example in ITtTI . 

Hence, we now present the fundamental mathematical description of the under- 
lying Monte Carlo algorithm that is developed for all the risk and insurance ap- 
plications above for specific models. More generally, MCMC algorithms are based 
on sampling a Markov chain with invariant measure 77. In this context, the lim- 
iting measure 77 is often called the target measure. It is not difficult to construct 
these random processes. For instance, let us assume that the law of X is reversible 
w.rt. some Markov transition K{x,dy). In this case, starting from the set A, we 
sample a sequence of random states using the Markov proposal K, rejecting se- 
quentially all the states falling outside the set A. The algorithm is well defined as 
soon as K{x,A) = K{Ia){x) > 0, and the resulting Markov chain X„ coincides with 
the Metropolis-Hasting algorithm with probability transition given by the following 
formulae 

M{x,dy) ■.= K{x,dy) U{y)+(^l-J K{x,dz) 5,{dy)- 

It is not difficult to check that r\ is an invariant measure of the chain with transi- 
tion M, that is we have that 

{r\M){dy) j r]{dx) M{x,dy) = ri{dy). 

The exact acceptance-rejection method discussed above and in (fT2] i corresponds 
to the special case 

Kix,dy)=F{X edy) 

In more general situations, the proposal transition K{x,dy) amounts of moving 
randomly around the starting point x. The individual (sometimes also called the 
walker) makes a number of tentative steps until it succeeds to enter into the desired 
set A. In general, the random state at that (random) hitting time of A is not distributed 
according to 77 . Roughly speaking, when the proposal transition K is based on local 
moves, the individual tends to hit the set A near the boundary of A. To be more 
precise, starting from an initial state Xq—xgW^—A the hitting time 

T :=inf{« >Q : X„eA} 

is a geometric random variable with distribution 

P(r=n \Xo=x) = {l~K{x,A))"'^ K{x,A) 

and we have 

E(/(Xr) I Xo=x)^KAif){x) := /:(/U)(x)//:(U)(x). 

When the chain enters in A, it remains for all times confined to the set A. In 
addition, under some weak regularity conditions on the Markov transition K, the 
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target measure rj is approximated by the occupation measures of the states; that is, 
we have the following asymptotic convergence result 

-TT L 5z„^„t~^? and ¥{X„edy\Xo^x):=M"ix,dy)^„^^ri{dy). 

" + ^ 0<p<n 

(16) 

In the above display, M" {x, dy) stands for the n compositions of the integral operator 
M defined by the induction formulae 

M"{x,dy)= J M"-\x,dz)M{z,dy) = j M{x,dz)M"-\z4y) 

with the convention M*'(x,c/3') = 5x{dy), for n = 0. 

It is of course out of the scope of this article to prove the ergodic theorem stated 
in the l.h.s. of ( fTSI l. We end this section with a simple proof of the r.h.s. assertion. 
Firstly, we observe that 

M"{f){x) = E(/(X„) lr<„ I Xo -x)+E(/(X„) It>„ \ Xq=x) 
= {\~K{xA)f-'K{xA)KA{M''-'{f)){x) 

l<k<n 

\k-l 



k>n 



On the other hand, we have 



lA{x)M{x4y) = \A{x)M{x4y)\A(j) M{\a){x) = \a{x) 

=^ \A{x)M{x4y) = \A{x)MA{x4y) 

with the Markov transitions 

M{\a){x) 

This clearly imphes that 

lAM"\f)^\AM'^{f)^KAM'"=KAM';i and r]M = r]M'l 
from which we find that 

KAM-{f){x)^r]{f)^ I KA[x,dy)r]{dy') [M;r(/)(y) -M^ (/)(/)] . 
This implies that 

sup||/^^M'"(x,.)-7?IL<^(m;;'):= sup \\M2{y..)-M2{y',.)\l, 
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In the above display, jj/ii — /X2||fv stands for the total variation distance between two 
probability measures jXi , and ji2 defined by 

||Ml-M2||n. = SUp{[jLii-jLi2](/) : osc(/)<l}. 

We consider the following mixing condition 

{Ha) There exists some probability measure v on R^', and some e]0, 1] such 
that 

Vx G A Ma (x, dy) > Ea v{dy) 
This condition is clearly met as soon as 

VxeA K{x,dy)lA{y)>eAvidy) 

For instance, when K{x, dy) has a density k{x,y) w.r.t. the Lebesgue measure X (dy) 
on R'', the condition {Ha) is met as soon as k{y) := mfx£Ak{x,y) is s.t. X{k) > 0. In 
this case, (Ha) is met with Ea = X{k) and v{dy) °^ k{y) X{dy). 
We also recall that 

(//^)=^j3(MX)<(l-e^r. (17) 

Next, we provide a short proof ( fTTI i. Under the rh.s. condition, the following Markov 
transition 

M',{x,dy)---'^^^^'^^^^^^^^ 
I-Ea 

is well defined, and we have 

MA{f){y) ^MA{f ){y') = (1 - £a) {MMy) -M'^{f){y')) ^ ^a) < (1 - e). 
Iterating the argument, we readily prove that for any m > 1 
sup \\M'^{y,.)~M2{y',.)\l<{l-eAr. 

Using the decomposition 

M"(/)(x)-77(/)= {\-K{x,A)f-'K{x,A) (KA{M"->^{mx)-n{f)] 

\<k<n 

+ (.m-riif)) I (1 -KixA))'-'Kix,A) 

k>n 

we prove that 

\\M^{x,.)-i]\\,,.< ^ {\-K{x,A)f-'K{x,A){l-£A)''-' 

\<k<n 

+ {l-K{xA)f-'K{x,A). 

k>n 
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After some elementary computations, we conclude that lim„-|^„ \\M" (;c, . ) — 77 = 0. 



3.4 Sequential Monte Carlo for Risk and Insurance 

The application of Sequential Monte Carlo (SMC) methods in the study of important 
problems in risk and insurance modeling is still relatively under developed, hence 
the motivation for this article. In the context of risk modeling see the example in 
||2l and the references therein for more discussion. We start this section with a mo- 
tivating class of algorithms targeting rare-event simulation via the restriction of a 
target measure to a contracting, increasingly rare set, such as a tail event. We detail 
the formal mathematical description of such a particle system stochastic algorithm, 
which is often omitted in the application papers utilizing such methods and some 
times to the detriment of both understanding and validity of results obtained. We 
believe it is important to present the principled understanding of this class of nu- 
merical integration regimes to clearly understand their utility in solving inferential 
problems in risk and insurance settings. 

However, we note that whilst this is a good motivation for these methods in risk 
and insurance, it is by no means the only method of SMC in risk and insurance 
applications as we illustrate in the class of algorithm developed in Section\^ 

Sequential Monte Carlo methods are acceptance-rejection techniques equipped 
with a recycling mechanism that allows to sample gradually a population of individ- 
uals w.rt. a sequence of probabilities with increasing complexity. We illustrate this 
methodology in the situation discussed above. Let us choose a decreasing sequence 
of subsets (Ap)o<p<„ joining Ao = R'' to the desired lower subset A„ = A: 

Ao =K'' cAi C A2 C ... CA„_i CA„ =A. 

Now, let's try to sample sequentially random copies of the random variable X w.rt 
the conditioning events {X G Ap}, with p < n. To get one step further, we let r\p be 
the sequence of measures 

7]p{dy) ■.= f{Xedx\Xe Ap) with p<n. 

By construction, (J7p)o<p<„ is a decreasing sequence of measures w.r.t. the abso- 
lutely continuous partial order relation pL between probability measures Q; that 
is, we have that 

nn < < . . . < 772 < 171 < rjo = Law(X). 



' we recall that /i v as soon as v(A) = /i(A) = 0, for all measurable subset A C M'' 
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3.4.1 Sequential Markov chain Monte Carlo methods 
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In this connection, we further assume that we have a dedicated MCMC style prob- 
ability transitions Mp with invariant measure r]p = rjpMp, for any p < n. We start 
running a sequence of random states {Xp)o<p<i,^ with transitions Mi, and initial 
condition tjo- For a sufficiently large time horizon «i, both the occupation measure 
■^T.i<p<n, ^Xp and the law of the terminal state Law(X„j) — rjoM"^ := Tti approx- 
imate the target measure rji. We also notice that the chain {Xp)p^<p<ni is confined 
to the set Ai as soon as one of the random states Xp^ e Ai hits the set Ai for some 
Pi < «!■ 

In the second step, starting from X„^ we run a sequence of random states 
(X„,+p)o<p<„, with transitions M2 (and initial condition n^). For a sufficiently large 
time horizon 112, both the occupation measure ^ Li<;j<ni ^x„^+p and the law of the 
terminal state Law(X„j+„, ) = K\M'^ approximate the target measure 772. As before, 
the chain (X„,+p)p2<p<n2 is confined to the set A2 as soon as one the random states 
Z„j+p2 ^2 hits the set A2 for some p2 < «2- 



JWfl A<2 iV/'S M?3 

1 



Tjo '- > rjoM"^ := tti ^ ttiM"^ = 712 > K2M"^ = Kt, 



(18) 



3.4.2 Feynman-Kac models 

Our next objective is to better understand the evolution of the flow of measures r\p, 
from the origin p = up to the final time p = n. Firstly, it is readily checked that 



'(X e dx I XeAp+,) = , 1^^. ^ W <^dx\X€Ap) 



1 

'{X€Ap~\X€Ap) 
and 

P (X e Ap+i I X e Ap) = J u^^, (x) F{xedx\xeAp). 

Therefore, in a more synthetic way, if we set Gp{x) = 1^^+, (x), then we have that 

rip+i = 1'G,XVp) 
with the Boltzmann-Gibbs transformation defined by : 



rip{dx) ^G„irip){dx) ^^-^ Gp(x) 77p(^^x). 

The next formula provides an interpretation of the Boltzmann-Gibbs transformation 
in terms of a nonlinear Markov transport equation 
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with the Markov transition S'p jjp defined below 

Sp,r,^{x,dy) = G„ix) 5.,{dy) + (1 - G„{x)) ^G„{rip)idy). 

In summary, we have shown that (T7p)o<p<« satisfies the following evolution equa- 
tion 

^OriQ '^l-Tl '^2.Tl2 "^3.^3 

i?o > rii > m > ri3 — ^ ^74 • ■ ■ 

In other words, rjp = Law(Xp) can be interpreted as the laws of the random states 
of a Markov chain X* with transitions Sp^ri,,', that is, we have that 

SO.TJO 59.7), ■^3.7)3 

Xq* > Xj* > ^ X3* ^ — > . . . 

The Markov chain X* can be interpreted as the optimal sequential acceptance- 
rejection scheme along the non increasing sequence of subsets Ap, in the sense that 

^p £ ^p+i ^ ^p+i ~ ^p 
X* e Ap-Ap+i =^ = 

where ^^^1 stand for an independent random sample with distribution 77^+1 = 
^Gpi^lp)- When the sample X* is not in the desired subset A^+i, it jumps instantly 
to a new state X**^^ randomly chosen with the desired distribution Tjp+i = fc^ (rjp). 
Next we provide a brief discussion on the optimality property of this Markov chain 
model. We recall that 

ll'?,,+ l-'?pl|rv 

= sup{[77p+i -77,J(/) : osc(/) < 1} 

= mf{¥{Xp^X,,+ i) : {Xp,Xp+i) s.t. Law(Xp) = Tjp and Law(Xp+i) = 77p+i } 

In the above display osc(/) = sup^--|,(|/(x) — f{y)\) stands for the oscillation of a 
given function / on K''. In this situation, it is instructive to observe that 

\\rip+i-rip\\n. = r{x;^x;+i). (19) 

In other words, the chain X* with Markov transitions S'p.Tjp realizes the optimal 
coupling between the sequence of distributions rjp. From the above discussion, we 
clearly have that 

P ^ x;) = Tjp(Ap - Ap+i) = 77p(l - Gp) = 1 - Tjp(Gp) 

On the other hand, we have 
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rip+i{f)-rip{f)^rip{S,,^^{f)-f) 

= '?p([l-Gp][/-fG, (%)(/)]) 

Choosing / = 1 — Gp, so that 

^G„ilp)if) = '^-1'G,irip){Gp)^0 

and 

rj, ([1 - Gp] [/ ~ ^G„ (r7p)(/)] ) = rip{[l- G„f) = 1 - 77,, (G,) 
This ends the proof of the optimal coupling formulae iT% . Next, we observe that 

l-rip{Gp) = l-rio{Ap+i)/rio{Ap) (with 77o = Law(X)) 

from which we conclude that 

i?o (Ap) > 770 (Ap+i) >(!-£) T?o (Ap) =^ P = x;) > 1 - e (20) 

As the reader may have certainly noticed, the Markov chain has very poor stabil- 
ity properties, in the sense that the distributions rjp strongly depends on the initial 
distribution 770 • More precisely, rjp coincides with the restriction of 7jo to the subset 
Ap-, more formally, we have that 

rip{dx) = fG„_, (7jo) = ^^^^^ ^ lA„ix) rioidx) 

The sequential Monte Carlo methodology is based on combining the MCMC 
methodology presented ( fTSb with the sequential acceptance-rejection technique dis- 
cussed above. To describe with some precision this method, we let Mp be an MCMC 
transition with invariant measure rjp ~ rjpMp. In this case, we have the evolution 
equation 

rip+\ = r]p+iMp+i = WG^{r]p)Mp+i := ^p+i(7jp) 

Notice that <Pp+i maps the set of probability measures 77 s.t. r\{Gp) > into the 
set of probability measures, and it is the composition of an updating transformation 
and a Markov transport equation w.rt. Mp+i ; that is, we have that 

Ip > Ip ■■= 'f'G„(j?p) > ^pMp+i = %+iirip) 

The solution of this equation is given by the Feynman-Kac measures defined for any 
measurable function / on IR'^ by the following formulae 

'?/>(/) = 7p(/)/7p(1) with 7p(/)=e(/(Xp) n G,{X,)]. (21) 

V 0<i/<p / 

To prove this claim, we use the Markov property to check that 
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rp+iif)=ElMp+iif){Xp)G„iXp) n GM]=Yp{GpMp+iif)). 

\ 0<q<p J 

This clearly implies that 

_ yp(G,M,+i(/))/yp(l) _ lpiGpMp+,{f)) 

^""^^^^ ~ 7piGp)/rpii) - npiGp) ~ ^-'^(^^)^''+'(^)- 

We already mention that the unnormalized measures 7,, can be expressed in terms 
of the flow of measures (T7p)o<p<n with the following multiplicative formulae 

7p(/)-i?p(/)x n "I'liG,) (22) 

0<c/<p 

This result is a direct consequence of the following observation 

7p(1)=e(g„_i(Xp_i) n =7p-i(Gp-i)-J7p-i(Gp-i)rp-i(l)- 

V 0<9<p-l / 

It is readily checked that the measures 77„ are the n-th time marginals of the 
Feynman-Kac measures on path space defined by the following formulae 



1 



n Gp{Xp) \ d^„ (23) 

0<p<n I 



with some normalizing constants = 7, ( 1 ) and the reference measures 

P„=Law(Xo,...,X„). 

This class of path space measures goes beyond the MCMC model discussed above. 
These measures represent the distribution of the trajectories of a reference Markov 
process, weighted by a collection of potential functions. These functional models are 
natural mathematical extensions of the traditional change of probability measures, 
commonly used in importance sampling. 

From a pure probabilistic viewpoint, these measures can be interpreted as the 
conditional distribution of a given Markov chain w.r.t. to a sequence of events. For 
instance, if we take G„ = Ia„ indicator potential functions of some measurable sub- 
sets A„ G E„, then we readily check that 

Q„ = Law((Xo, . . . ,X„) I VO < /9 < n Xp e Ap) and JJ, = P (VO < p < « Xp e Ap) 

In filtering settings, if we take G„{x„) = p{y„\x„) the likelihood function associated 
with the observation Y„ = y„ of the random signal state X„, then we have 

Q„=Lav^{{Xo,...,X„)\yO<p<n Yp^yp) and X ^ p{yo, ■ ■ ■ ,yn-i) 
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For a more thorough discussion on the appHcation domains of these Feynman-Kac 
models, we refer the reader to the books 1261 |23l l24l [3T1 . 

3.5 Nonlinear distribution flows 

The central idea behind Feynman-Kac particle samplers is to observe that any evo- 
lution equation of probability measures 

on some measurable state spaces £„ can be interpreted as the law 

T],, = Law (Xn) 

of a Markov chain X„ with initial distribution Tjo and Markov transitions 

G dx„ I =x„_i) — Ki,^ri„_[{xn-iTdx„). 

The Markov transitions K„ r\„_^ are chosen so that 

V« > 1 7?„-l^r„,j7„_i = <?„(t7„-i). 

The Markov chain Xn incorporate free evolution moves according to M„, with se- 
quential updates of the measures 7],,, so that the law of the random states X„ coincide 
with the desired distributions r\„, at every time step. In this interpretation, this chain 
can be interpreted as a perfect sequential sampler of the sequence of measures T],,. 

The choice of the transitions /C„+i jj„ is not unique. For instance, for the Feynman- 
Kac models on E„ = W' discussed above, if we take 

Kn+\,T],X^,dy) := [5„,,,„M„+i] {x,dy) or K^+i^ri.X^idy) := (tj„) (dy) 

we readily check that 

rinK„+lji„ = <P„+l ()?„) = 'FG,Xnn)M„+l = rinSn.n„Mri+l- 

We also mention that the law of the random trajectories {Xq, . . . ,X„) are given by 
the so-called McKean measures 

F„{d{xo,...,Xn)) = rio{dxo) Ki,j^^{xo,dxi) . . .K„^ri„_i{^n-i,dx„). 

We further assume that the Markov transitions M„(x„^i,dx„) are absolutely con- 
tinuous with respect to some reference measure v„ and we set 

Qn{x„_i,dx„) := G„^i{x„_i)M„{x„-i,dx„) = //„(x„_i,jc„) V„{dx„). 

In this situation, we have the following time reversal formulae 
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Q„(^f (xo, . . . ,x„)) — r}n{dxn) M„^,]^_j {x„, dxn-i)...Mi^^^{xudxo) (24) 
with the Markov transitions 

^„ / J \ ^n-l{dx„^i) Hn{Xn-l,Xn) 

We prove this backward formula using the fact that 

ri„{dx„) = fG„_|(Tj„_i)M„(c/x„) = ^"-1 {Hn{-,x„)) ^^^^^^^ 

from which we find that 

r]n[dXn) Mn.n„_i{Xn,dXn-\) = JT:. T ^n-\{dXn-\) Qn{Xn-\,dx„) . 

Iterating this process, we prove (|24] |. 



3.5.1 Interacting particle methods 

This section is concerned with particle approximations of the Feynman-Kac model 
(I2TI 1 and ( l23T l. We also present a series of exponential concentration inequalities 
that allows to estimate the deviation of the particle estimates around their limiting 
values. 

In the further development of this section /„ stands for some function s.t. ||/„ j| < 
1, and (ci,C2) represent two constants related to the bias and the variance of the 
particle approximation scheme, and c stands for some universal constant. The values 
of these constants may vary from line to line but they don't depend on the time 
horizon. Last, but not least, we assume that the Feynman-Kac model satisfies some 
strong stability properties. For a more detailed description of the stability properties, 
and the description of the quantities (c,ci ,C2) in terms of the Feynman-Kac model 
(|2T), we refer the reader to the books 12311241 . 

We approximate the transitions 

Xn -^Xn+\ ^ Kf,+ij\n{Xn,dXn+\) 

by running a Markov chain ^„ ~ (4,' ,...,t,^)<^ that approximate the distribution 
7],, when N^°° 

]^ L % ■= -^m^ nn. 

" l<i<N 

A natural choice of particle transitions is to take at every time step sequence of 
conditionally independent particles 
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For the Feynman-Kac models discussed above, we can chose the transitions Ar„+i .r]„ = 
Sn,r]„M„^ 1 • In this context, the evolution of the particle algorithm is decomposed into 
two steps. 



4,' 



+1 



+1 



During the first step, every particle 4i evolves to a new particle randomly chosen 
with the distribution 

S^.{^;„dx) G„(4;) 5^,{dx) + (1 - G„(^i)) -ffe^lOl^^) 
with the updated measures 



This transition can be interpreted as an acceptance-rejection scheme with a recycling 
mechanism. In the second step, the selected particles 4', evolve randomly according 
to the Markov transitions M„+i. In other words, for any 1 < i < A^, we sample a 

random state ^^^^j with distribution M„+i (^^^,dx^. 

Using the concentration analysis of mean field particle models developed in ll72ll . 
the following exponential estimate was proved in 124]. For any x > 0, « > 0, and any 
population size N >l, the probability of the event 



C2 



is greater than 1 — e In addition, for any x = (j^!)i<,<c/ and any (—°°,x] = 
nf^j {-oo^Xi] cells in E„ = R'', we let 

F„{x) = ri„ and (x) = 77,^ . 

For any y>Q,n>Q, and any population size > 1 , the probability of the following 
event 

Vn \\F„''~F„\\<c^d{y+l) 

is greater than 1 — e^-'. 

If we interpret the mutation-selection particle algorithm as a birth and death 
branching process, then we can trace back in time the whole ancestral line Q = 
i^'p n)o<p<n of the individual at the n-th generation. 
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The random state „ represents the ancestor of the individual at the level p, with 
< p <n, and 1 < / < A^. It is more or less well known that ^„ coincides with the 
particle approximation of the Feynman-Kac model defined as in ( l2Tl i by replacing 
X„ by the historical process (^p)o<p<n- This interpretation provides an alternative 
particle approximation scheme of the measures (|23] |. that is we have that 



More precisely, we proved in 11241 the following exponential concentration estimate. 
For any test function fn on path space s.t. ||fn|| < 1, for any y > 0, « > 0, and any 
N > I, the probability of the event 



K - Q„] (/) < ^ (1 +. + V^)+C2^J ^ 

is greater than 1 — e^''. 

Further details on these genealogical tree models can be found in 1231 l24l l73l . 
Mimicking formulae (|22] ) and ( l24b . we define an unbiased particle estimate "/^ of 
the unnormahzed measures 7„ and a particle backward measures by setting 

i^(/)='7:(/)x n <(G.) 

0<c/<n 



and 



{d{xo, . . . ,x„)) = 7]^ [dxn) ^n.-^N ^ (x„,c/x„_i) . . .Mj (xi ,c/xo) 



We end this section with a couple of exponential concentration estimates proved 
in im . For any jc > 0, « > 0, > 1, and any e e { + 1,-1}, the probability of the 
event 

^log^<^(l+x+VJ)+^V^ 
n ^ 7„(1) ^ ' ^ff^ 

is greater than 1 — e^"^. In addition, for any normalized additive functional fn(.xo, . . . ,x„ 
;riT^o<p<n//)(^p) with ll/pll < 1, forx > 0, « > 0, and any population size > 1, 
the probability of the event 



m - (f") < 1 (1 + (x + + C2 y^l^ 

is greater than 1 — e"'^. 
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4 Illustration of Interacting Particle Solutions for Risk and 
Insurance Capital Estimation 

In this section we detail a special sub-set of algorithms from within the stochastic 
particle integration methods that was specifically developed to solve problems for 
risk and insurance in |[36l and discussed in comparison to specific FFT methods in 
||63]| . The class of recursive solutions developed is very general and applicable to a 
wide range of insurance and risk settings. We provide a novel result in this illustra- 
tion which extends the framework originally presented in 1361 through consideration 
of a higher-order recursion based on the widely utilized Panjer recursion (avoiding 
the need to perform discretization of the severity distribution). 

It is important to understand that whilst this method is focusing on asymptotic re- 
sults for tail functionals of a compound processes model under the LDA framework 
and can be re-engineered into a class of rare-event type algorithm as discussed in 
Section 13.41 the presentation of this example that we adopt will be more general 
and therefore more applicable to a wide class of insurance and risk problems that 
include but is not limited to those used to motivate the use of stochastic particle 
solutions in the introduction. 

Finally, we point out that this method to be discussed is one of many methods that 
could be developed for such an insurance problem, though it has proven to be highly 
efficient even in the simplest form of the sampler as discussed below. In addition we 
note that as discussed, whilst the asymptotic approximation methods that were pre- 
sented were developed to tackle the serious statistical and computational challenges 
posed by accurate estimation of tail quantiles and expectations for heavy-tailed LDA 
models. Noting that their elegance lies in the fact they allow one to bypass the seri- 
ous computational challenges for estimation of risk measures for such heavy-tailed 
annual loss distributions under traditional integration methods, Fourier methods, re- 
cursions (Panjer) or basic Monte Carlo approaches, they do have associated issues, 
see discussions in ||56l and the arguments presented in the introduction of this ar- 
ticle. As discussed the properties of such SLA estimates is still an active subject 
of study with regard to approximation error, unbiased quantile function estimation, 
asymptotic rates of convergence, sensitivity to parameter estimation and model mis- 
specification. Hence, one often requires calculation of VaR, ES and Spectral Risk 
Measures which do not utilize such asymptotic properties but instead are based on 
stochastic integration methods. Here we consider such an example. 



4.1 Recursions for Loss Distributions: Panjer and Beyond 

The framework proposed in ll36l and ||74| for developing a recursive numerical so- 
lution to estimation of such risk measures through estimation of the density of the 
compound process. In particular we briefly summarize an approach to transform 
the standard actuarial solution known as the Panjer recursion llTSl to a sequence of 
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expectations. Recursions for the evaluation of single risk process distributions are 
ubiquitous in risk and insurance modeling, see the book length review in discrete 
loss distributional settings in fTS^. In this paper we consider an advanced develop- 
ment that avoids the need to discretize the severity distribution to perform inference 
and instead develops a stochastic particle integration based solution, hence provid- 
ing a fitting example for such methods in risk and insurance. 

Recursions for evaluating Fz{x) and risk measures for Subexponential LDA 
models are traditionally in risk modeling based around in the most fundamental 
class of models, the Panjer class Panjer Recursion Class of frequency distribution 
relationships defined by, 

Pr,= ( « + - ) Pn-l (25) 



n 

with members Poisson {a = 0,b = l,po^ e^^). Binomial (a = jf^i^^ ^'("i -g/ ' ~ 

(l—q)m) and Negative Binomial (a = j^,b^ TTIT'/'o ~ {l+b) — r).ln addition, 
one can can derive the generalized (higher order Panjer) recursion for an extended 
class of frequency distributions given by the generaUzed Poisson distribution (GPD) 
via 




V« = 0,l,2,... 

if « > m, when < 0. 



with X >0 and max(— 1 , X /m) < < I and m > 4 is the largest positive integer s.t. 
X + dm > when Q is negative, where the GPD is Poisson for 0=0; over-dispersed 
> and under-dispersed < 0. 

Considering these classes one can derive closed form recursions for the annual 
loss LDA compound process distribution according to the Panjer recursion ifTSl , 

/f (x) ^infx{x)+ (« + y) ^ - (26) 



or the generalized higher order Panjer recursion 1771 , 

A 



fy {x) = p 1 (A , )fx (x) + j-^ J^(^e + Xl)fx (y) fr {x - y) dy 



(27) 



To understand where these recursions are obtained from, consider the convolution 
identity for an i.i.d. partial sum S„+\ = X\ + . . . +X„+i with density 

f f{z)r{x~z)dz, V«= 1,2,3,.... (28) 
Jo 

Substitute the conditional of X\ when 5„+i = x. 



(29) 
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into the average given 5„+i = x to get 

E[X,\X, + ---+X„+,=.x]= ' dr. (30) 



Then observe that with i.i.d. losses one also gets 

I n+l 

E [Xi \Xi + --- +X„+i = x] = —— ^ E [Xi\Xi + • • • +X„+i - x] 

n + 1 ,t1 

= [Xi + -- ■+X„+i\Xi + ■■■ +X„+i = x] = ^ 

n+l n+l 

Equating these conditional mean expressions and rearranging gives 

n + l ' X Jo ' 



(31) 



(32) 



Now utilize the Panjer class of frequency distributions satisfying for some a and b, 

FT{N = n)=p„=(^a + ^^Pn-u (33) 
Upon substitution and some elementary algebra one obtains the Panjer recursion. 



4.2 Stochastic Particle Integration Methods as Solutions to Panjer 
Recursions 

Solving these recursions is then typically performed to obtain estimates of Fz^^ (x). 
In the case of sub-exponential severity models this can be costly if the standard 
actuarial approach is adopted which involves discretization of the severity distribu- 
tion. This renders the integral equation in a form for certain classes of frequency 
and severity models that admit the De Pril transform solution or approximations of 
such a transform solution, see discussion in 1761 . Such approaches are then prone 
to a trade-off between discretization errors and computational efficiency, though the 
errors are deterministic. The second approach one may consider is to avoid the dis- 
cretization error in favor of a Monte Carlo solution. If this can be done efficiently 
with a clear understanding of the associated Monte Carlo errors, this may be con- 
sidered as a preferred alternative. Such an approach was adopted under a path space 
based Importance Sampling solution in [361 . 

It was noted in 1361 that since the Panjer recursions could be re-expressed as 
linear Volterra integral equations of the second kind via the mapping 
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xi =x~y, g{x) = pifx{x), f{xi)^fz{xi), and 



X—Xl 



k{x,xi) = a + b ] fx{x — xi) . 



(34) 



where the kernel k : E x E i-^M. and function g : E i-^M. are known and / : £ i— >■ R 
is unknown. Furthermore, if one defines k''''{x,y) = 1, k^ (^jj) — K^^y) ™d 

k"{x,y)^ j k{x,y)l^-\z,y)dz 

and these kernels satisfy that 

^ / {xQ,x„)g{x„)\dx„<oo 

then one can identify the resolvent kernel and Neumann series through iterative 
expansion of the recursion to obtain for a sequence of domains E\;„ 

f (xo) = gixo) + ■■■/ g{xn)Y[k{xt_i,xi)dxi;„, 

„^qJQ jo 

see further details in ||36l . Under this formulation it was shown in [36:] how to ad- 
dress two problems: estimation of the annual loss density over a set A and estimation 
of the annual loss density point wise. These are both directly relevant to obtaining 
estimates of the risk measures specified. 

To proceed one converts the Neumann series into a sequence of expectations with 
respect to an importance sampling distribution given by associating the following 
elements 

n 

fo{xo) =g{xQ), andf„{xo:„)=g{xn)Y[k{xi^i,xi) 

1=1 

■' ■ f (xo) = fo (xo) + y ■■■ fn{xO:n)dxo:„. 

,^iJo Jo 

in order to frame this problem as an expectation with respect to a sequence of dis- 
tributions {7r(«,jci:„)}„>o: 

w X /oW^m\ , V f fn{x,Xi;„) 

■^W = :;;77TV^ + L / , , —f rn{n,xi;n)dxi:„ 



7r(n,xi:„) 



with the sets Ai:„(,yo) = {{xi, . . . ,x„) : xq > xi > ■ ■ ■ > x,,} playing an analogous 
role to the sequence of level sets described in the methodology section of the pa- 
per. 
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Hence we note there are clearly now two path-space based particle solutions 
available, those that consider estimating f{x) point-wise which involves an impor- 
tance sampling solution on the path-space defined by 

Ur=o{«}x^l:nW- 

Note, efficient variance reduction (per sample basis and per unit of computation ba- 
sis) involves instead estimating (/(x) — /o (xq)) by importance sampling on smaller 
space 

U"=l{«}xAi:„(x). 

The other alternative involves characterizing f{x) over some interval by obtaining 
samples from its restriction to that interval [x„,x;,], via importance sampling on a 
slightly larger space 

U {«} xAi:„([x„,Xi]). 

One can now consider constructing importance sampling based solutions to this 
sequence of expectations as detailed in 1361 [Algorithm 1, p. 9] and [Algorithm 2, 
p. 12] and |[74l [Algorithm 2. 1.1] and for which we provide the relevant pseudo code 
in Algorithm 2. This is summarized according to the path-space based Sequential 
Importance Sampling (SIS) based approximation to the annual loss distribution, for 
a Markov chain with initial distribution/density jU (x) > on £ and transition kernel 
M{x,y) > if ^(x,^) 7^ and M has absorbing state d such that M{x,d) = Pj 
for any x e £, by the following steps in Algorithm 2. 

Algorithm 2: Stochastic Particle Methods on Path-Space for Panjer Recursions 

1. Generate independent Markov chain paths absorption 

«(') + ! 

2. Evaluate the importance weights for each particle on the path space by. 



1 



W(Aj = <C V ^/ (35) 



0:«(' 



4.3 Stochastic Particle Solutions to Risk Measure Estimation 

Then if /i ( ^q'' ) = ^ ( ^d'' ) ? ^' £ {!:••• then the empirical measure at a point 



xo is given by 



/z(-o) = ^Ew(xo,X«„ 
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or over an interval, 



/z(-o) = ^Ew,(x»,,)5(xo-4 

(=1 



forms an unbiased Monte Carlo approximation of the expectation of fz{z) for any 
set A given by E jj^ f{xQ)dxQ ~ J^f{xo)dxo- Furthermore, detailed discussions 
on the optimal choice with respect to minimizing the variance of the importance 
weights is developed in ll36l and ll74l . 

Having obtained this particle based approximation, this weighted dirac measure 
can then be utilized to estimate any of the required risk measures such as VaR, ES 
and SRM for any desired level a. This can be performed in two ways, depending on 
whether the particle solution is obtained for the evaluation of the recursions point 
wise over a fixed grid or alternatively over an interval, which could be increasing 
in size. In the case of an interval, or contracting set, one considers perhaps a set 

of interest to be A = [0,Xmax] such that x„„„ >> ^1 — ^p^^, and then utilize 
this to construct an unbiased particle approximation of the distribution of the annual 
loss up to any level a E (0, 1). Clearly, this could be obtained from growing a set 
Ai = [OjXi] C A2 C • • • C A = [OjXmax] recursively, as discussed in previous sections 
where such as set would result in the tail F being restricted to a smaller and smaller 
set of the complement rare event. 

If the partition is made point wise over a linear or a non-linear spacing [Q,z] = 
Um=i [('"^ l)A,mA) and the distribution evaluated point-wise, this leads to an 
estimation of 

Fz{z) = E A/(mA) « - E E AW (mA,^;;;;;:^ . (36) 

m=0 777=0 7= 1 

Alternatively, if the estimation is performed over an interval, this construction of 
the empirical measure over interval A(xmax) ~ [0,Xmax] such that for any z < Xmax 
results in. 



^z(^) = ^fw «i„) e [0,z]) ^^t» /V^W^^- (37) 



Practical advice: consider a range for the support [0, Xmax] s.t. x,j,ax 

From these unbiased particle approximations of the annual loss density and dis- 
tribution, the evaluation of the risk measures for VaR, ES and SRM follows triv- 
ially. To begin with the VaR or indeed any quantile, it would be useful to have a 
reconstruction of the quantile function of the annual loss LDA model distribution 
function. The quantile function is defined by 

Qip) = F^ip) = inf {-^ eR: p< Fzix)} . (38) 



Title Suppressed Due to Excessive Length 37 

Estimation of the Quantile function can now be obtained by approximation utilizing 
the empirical measure either based on a random set of particle locations or a discrete 
deterministic grid as follows: 

Deterministic Grid Solution: Given partition [0,Xmax] = Um=i [{m — 1) A,mA) 
for some step A s.t. 




Guide - set x,nax >> (l - ipf ) ■ 

Interval Solution: Construct the empirical measure over A(oo) = [0,°°) s.t. 

(40) 

XjQj represents the order statistics for the particles. 

From these distributional estimates for density, distribution and quantile function, 
one can then obtain a particle based weighted dirac measure can estimate the risk 
measures for any a e (0, 1) by: 

Value-at-Risk (VaR): directly obtained using the estimated quantile function! 
Expected Shortfall (ES): estimated via 

Spectral Risk (SRM): the SRM for a weight function ^ : [0, 1] R is given by 

SRM^('?') = ^E^(0):„(..'/'(/'')Ap, 

withp,=l,.,,iv(x«^,,,). 

For additional discussions and detailed examples and applications of this numerical 
approach to risk estimation can be found in 1361 and [?], and in financial asset 
pricing in (741 ■ To conclude this illustrative section, we extend the previous example 
based on the popular Poisson-Log Normal LDA risk process model to include the 
path-space particle solution discussed above. 

Example 2 (Poisson-Log Normal LDA Model (continued)). Consider the Poisson- 
Log Normal compound process detailed in Example 1 . We demonstrate results for 
standard Monte Carlo approach and compare results to the path-space particle so- 
lution discussed and detailed in Algorithm 2. The Monte Carlo solution generated 
N=50mil samples (hence effectively exact as Monte Carlo error is insignificant) and 
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a grid based solution was adopted for the particle solution with N=50k per grid 
point giving a total of NT=500k, with a grid width = 1 . The results for the estima- 
tion of the quantiles (rounded to integer) which includes the Basel II/III regulatory 
standard for Economic Capital and Regulatory Capital reporting are provided in the 
following table for two sets of parameter settings of A = 2, /i = 2 and a = 0.5 and 
X — 2, pL = 2 and a = 1 . Two sets of parameters were considered, one based on 
the observed poor finite quantile performance of the Single Loss Approximate and 
the other in the setting in which the SLA is known to perform more accurately. The 
Monte Carlo errors for the Standard Monte Carlo approach are insignificant due 
to the massive sample size, hence we treat this as the exact solution upon which 
we compare the particle solution and the single loss approximation. The particle 
solution presents a 95% confidence interval and the single loss approximation sim- 
ply reports the asymptotic approximation no explicit error can be calculated for the 
point estimated quantile (as discussed above). 



Quantile Level 


Standard Monte Carlo 


Particle Solution (Alg. 2) 


Single Loss Approximation 




£7 = 0.5 


£7= 1 


£7 = 0.5 


£7= 1 


(7 = 0.5 


£7= 1 


50% 


14 


16 


15 [14,16] 


16 [13,17] 


10 


14 


80% 


27 


39 


25 [26,28] 


41 [39,43] 


14 


26 


90% 


35 


57 


33 [31,35] 


55 [52,59] 


16 


38 


95% 


42 


77 


40 [38,43] 


74 [70,79] 


19 


52 


99% 


57 


129 


55 [54,56] 


123 [119,127] 


26 


97 


99.5% 


77 


234 


73 [68,79] 


227 [218,240] 


38 


198 


99.95% 


83 


276 


79 [73,91] 


270 [261,282] 


42 


240 



Table 1 Standard Monte Carlo Solution (exact) versus Particle Solution and First Order Single 
Loss Approximations. 



The results that we present in Table I are obtained on a linearly spaced grid of width 
1 . However, this can be changed to either include a non-linear spacing, placing more 
points around the mode and less points in the tails or as we detailed straightout eval- 
uation on an interval, avoiding the discretization of the grid. For the sake of compar- 
ison between the standard Monte Carlo and the importance sampling estimates, we 
histogram the standard Monte Carlo procedure samples using unit length bins. We 
can see two things from Table I, firstly as expected the particle based solution per- 
forms accurately under any parameter settings for a modest computational budget. 
When compared to the Single Loss Approximation, we see that there is two clear ad- 
vantages in having a complementary particle solution, since we obtain measures of 
uncertainty in the quantile point estimates, trivially. Secondly, we demonstrate that 
the Single Loss Approximations may not be as accurate as required for even these 
simple models at quantiles that may be of interest to assessment and are required for 
reporting of capital figures under financial regulation standards. 
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